#ifndef _SUPPRESSOR_H_
#define _SUPPRESSOR_H_

#include <fftw3.h>
#include <vector>
#include <complex>
#include <algorithm>
#include <functional>
#include <numeric>
#include <iostream>

#include "criterion.h"

template
<
    typename Compression
>
struct Suppressor
{
public:

    template
    <
        typename SignalSystem
    >
    SignalSystem operator()(SignalSystem& sys, double compression = Compression().value)
    {
        using namespace std;
        
        auto result = sys;
        cout << result.size();
        for(int i = 0; i < result.size(); ++i)
        {
            cout << "signal size" << result[i].size() << endl;
            
            auto spec = Suppressor::fft(result[i]);
        
            cout << "spec size" << spec.size() << endl;
            // suppress
            fill(
                spec.begin()
                + min(
                    (double)spec.size(),
                    spec.size() * compression / 2.0
                ),
                spec.begin()
                + min(
                    (double)spec.size(),
                    spec.size() * (1.0 - compression / 2.0)
                ),
                0.0
            );
            
            cout << " a = " << min(
                    (double)spec.size(),
                    spec.size() * compression / 2.0
                );
            cout << " b = " << min(
                    (double)spec.size(),
                    spec.size() * (1.0 - compression / 2.0)
                ) << endl;
                
            
            result[i] = Suppressor::ifft(spec);
            copy(result[i].begin(), result[i].end(), ostream_iterator<complex<double>>(cout, " "));
            cout << endl;
            copy(spec.begin(), spec.end(), ostream_iterator<complex<double>>(cout, " "));
            
            cout << endl;
            cout << "signal size" << result[i].size() << endl;
            
            cout << "spec size" << spec.size() << endl;
            
            cout << endl << endl;
        }
    
        return result;
    }

private:

    static std::vector<std::complex<double> > fft(const std::vector<std::complex<double> >& signal)
    {
        using namespace std;
        
        std::vector<complex<double> > result(signal.size());

        // Create plan
        fftw_plan plan = fftw_plan_dft_1d(
            min(signal.size(), result.size()),
            (fftw_complex*)(&signal[0]),
            (fftw_complex*)(&result[0]),
            FFTW_FORWARD,
            FFTW_ESTIMATE);

        // Forward DFT
        fftw_execute(plan);
    
        // Free memory
        fftw_destroy_plan(plan);
    
        return result;
    }

    static std::vector<std::complex<double> > ifft(const std::vector<std::complex<double> >& signal)
    {
        using namespace std;
        
        std::vector<complex<double> > result(signal.size());
        
        // Create plan
        fftw_plan plan = fftw_plan_dft_1d(
            min(signal.size(), result.size()),
            (fftw_complex*)(&signal[0]),
            (fftw_complex*)(&result[0]),
            FFTW_BACKWARD,
            FFTW_ESTIMATE);

        // Backward DFT
        fftw_execute(plan);
    
        // Free memory
        fftw_destroy_plan(plan);
    
        return result;
    }
    
    static double add(std::complex<double> x, double y)
    {
        auto n = abs(x);
        return n * n + y;
    }
};

#endif
